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ABSTRACT 


This report analyzes the use of the Cholesky decomposition tech- 
nique as applied to the feature selection and classification algorithms used 
in the analysis of remote sensing data (e. g. as in LARSYS), This tech- 
nique is approximately 30% faster in classification and a factor of 2-3 
faster in divergence, as compared with LARSYS. Also numerical stability 
and accuracy are slightly improved. Other methods necessary to deal 
with numerical stability problems are briefly discussed. 

It is, in our view, extremely important that the best numerical 
techniques be used in production calculations. The ar^ment that sub- 
optimal techniques have sufficed in the past is not valid if one considers 
that unexpected failures in the future may be extremely costly to rectify; 
since the use of the validated techniques discussed in this report are 
more reliable and efficient, it would seem wiser to proceed into further 
production calculations with the assurance that the systems and methods 
used rest on a more secure algorithmic foundation. 
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THE USE OF THE MODIFIED CHOLESKY DECOMPOSITION IN DIVERGENCE 


AND CLASSIFICATION CALCULATIONS 


1 . Introduction : 

This report analyzes the use of the Choi esky decomposition^^ tech- 
nique in the analysis of remotely sensed data, specifically in divergence 
calculations and in the evaluation of the maximum likelihood function; the 

latter occur in, respectively, the feature selection and classification 

( 2 ) 

techniques used, for example, in the LARSYS' system developed by the 
Laboratory for the Applications of Remote Sensing of Purdue University. 
Although LARSYS was primarily developed for research purposes, in- 

( 3 ) 

creasing use of the system and of derivative systems such as ERIPS for 
production processing emphasizes the need for efficient, accurate and stabl 
algorithms as the basis for design objectives of computational analysis. 

The organization of computation in certain segments of LARSYS and the use 
of subroutines such as MINV from the IBM Scientific Subroutine Package 

/ccn \ 

IbbKj do not lend credence that such design objectives have been met. 

The purpose of this report is to describe how one possible re-organization 
of the computation and the use of preferred techniques can improve the 
efficiency and accuracy of the system. 


*Numbers in superscripts refer to references 
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The focus of this report is on improved efficiency in terms of compu- 
tation time. Thus (Appendices A and B) the arithmetic precision used is 
identical with that used in LARSYS, so that a meaningful comparison of effi- 
ciency can be obtained. It will be shown that the algorithms proposed yield 
improvements in computational speed with no loss in accuracy or stability 
(in fact, slight improvements can be obtained in the latter). 

Improvements in accuracy and stability can be achieved by further re- 
finements in the techniques used. This will be the subject of a later re- 
port; however, in Section 6, we discuss where such improvements can be ex- 
pected by the use of higher precision and/or the use of such techniques as 
iterative refinement, scaling and equilibration. 

It is, in our view, extremely important that the best numerical techniques 
be used in production calculations. The argument that sub-optimal techniques 
have sufficed in the past is not valid if one considers that unexpected fail- 
ures in the future may be extremely costly to rectify; since the use of the 
validated techniques discussed in this report are both more reliable and 
efficient, it would seem wiser to proceed into future production calculations 
with the assurance that the systems and methods used rest on a more secure 
algorithmic foundation. 

2. ChOlesky Decomposition : 

Let K be real, nxn, symmetric positive-definite matrix. In the 
applications under consideration, K would be a covariance matrix. Then 
there is a unique, nxn, real, lower- triangular matrix, L, such that 
( Choi esky decompos i ti on ) 


K = LL* 


( 2 . 1 ) 
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where L* denotes the (conjugate) transpose of L. There is also a unique, 
real, lower triangular matrix, L, and a real, positive diagonal matrix, 

D, such that ( modified Cholesky decomposition ) 


K = LDL* 


( 2 . 2 ) 


where L has diagonal elements equal to unity. From (2.1) and (2.2) it 
can be seen that 


L = LD^ 


(2.3) 


/2 

where D ' is the diagonal matrix whose entries are the square roots of 
the corresponding elements of D. 

Either the Cholesky or modified Cholesky decompositions can be readily 
obtained from the following recurrence relationships^^ ^ (we use the 
notation K = (l^j^j)/ l = (^j^j)» L = (^^^j)* d = diag 
D = diag {^^j-) ; 


Cholesky 


*11 - 

^11 . , 
D-1 

,2 

A. . = 



J3 

\ gg- L 

S=1 

gs/ 

1. . = 


i.l. 

3-D 

\ ij 

IS js 


s=l 

i=j+l, j+2, . . . ,n 

and, of course, = 0 for j > i . 


j = 1, . . . ,n 


(2.4) 
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Modified Cholesky 



i=j+l, . . . ,n (2,5) 


where = 1 (i=l,...,n) and = 0 for j > i. 

For the applications under consideration, the modified Cholesky de- 
composition is more useful since it avoids the computation of square roots 
inherent in (2.4). It can easily be shown that, under the assumption that 
K is positive-definite, S'. > o (j=l,...,n). 

Once either decomposition is obtained, solutions of equations of the 

form 

Kx = b (2.6) 


may readily be obatined from the back and fonvard substitutions (we henceforth 
only consider the modified Cholesky decomposition): 


Yl = L 


(2.7) 


= D V 


( 2 . 8 ) 


X = L 


( 2 . 9 ) 



KX = LDL*X 



= Ly^ 
= b 


as desired. (2.7)- (2.9) may alternatively be written (using •“ to denote 
replacement as opposed to equality) to economize on storage: 



i=n-l,n-2, . . . , 1 


( 2 . 10 ) 


Note that in order to solve such systems there is no requirement to calcu- 
late K“\ only L and D which requires approximately 1/3 the amount of com- 
putation. 

This saving in itself is significant if one considers that the amount 
of time devoted to computing matrix inverses in connection with feature ex- 

O 

traction in LARSYS varies roughly as mnf , where m is the number of 

classes and n is the number of features under consideration - the corresponding 
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amount of time devoted to the actual divergence calculation varies as 
hm n , which is of the same order of magnitude for most problems considered. 
Thus reducing the first factor by a third can significantly effect the over- 
all computation time of itself. 

In the applications under consideration, we thus have m covariance 
matrices Kg (s=l ,... ,m) corresponding to the number of classes. The di- 
mensionality, n, of each corresponds to the number of channels. With 

obvious notation, we write 

/s/ 

K = Lj D L S— 1 , . • . , ITI 

s s s s 

where 


U<5>) . £3 - 


and^ . {df>} 


are calculated as in (2.5) 


3, Feature Selection : 

Feature selection, as implemented in LAR5YS, depends upon calculating 
a measure of inter-class divergence for multiple classes, requiring calcu- 
lations of the form 


D 


Di + 


(3.1) 
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where 


m-1 m 

^ ■ -.1 

3 = 1+1 


y tr[ (K.-K.) (K ^-K ^)1 (3.2) 

-I ^ ^ J J X J 


m-1 m 

D_ = y y (M -M)* (K."^+K ”^)) (M-M.) 
2 j=i:+l 1 3 1 D 13 


(3.3) 


where tr a denotes the trace of A (sum of its diagonal elements) and 

M (s=l,...,itO is the mean vector for the s^*^ class. We first simplify 
s 

(3.2) and (3.3). 

We note that we can wri te 


m -1 

m 

m -1 

m 

II 

y (trK.K."^) 
j=-I+.l ^ 

3=1 

V’ -1 

/ ^ (ti^ .K. ) - nm(m- 

1=3+1 

% 

m-1 

m . 

m 

i-1 

= I 

y (trK.K."^) 

+ / 

y (trK.K. - nm(in-i) 

ii'i 

3=1+1 


j£l ^ 3 

m 

m 




y tr(K.K."‘^) 

2 

- nm 


m 




II 

.1^ tr(K.-lKj) 

2 

- nm 


tr (AB) 

= tr(BA)) 




m 




K) - nm' 


(3.4) 
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v/here 

m 




K = ) K. 

ifel ^ 




= LDL* (say) 



Now 

1 — 1 />-/ — 1 —I 

trK. K = tr(L. D. L. LDL* 

D 3 3 3 




_ *1 

= tr (D . T . D T .*) 

3 3 3 



where 

T. = L “^L = (t. .) (say) 

J J -^J 


(3.5) 

Thus 





tr(K.'^K) = i y 

3 ' pfei pq p q 


(3.6) 


• m f in 




Hence D-| may economically be calculated from (3.4), (3.5) and (3.6). It 
should be noted that the calculation of the |Tjj in (3.5) each require n 
calculations of the form (2.10); however, since T j ,*1 j and T are all lower 
triangular, it is important to remark that much of the computation may be re- 
duced by observing that, in calculating the column of Tj , the index 
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n in (2.10) is actually replaced by n-q+1 (q=l,...,n). 

The calculation of may be similarily simplified. For, from (3.3), 
we may write 


m-1 m 

D, = ) ) {M-M ) * (K,"^+K.“^) (M-M.) 

2 ifrl j=t+l 1-313 13 


m-1 m 


= ) y (M-M )*K. ^(M -M ) 
i=l j=i+l ^ ^ 

m-1 m 
i=3+l 


(3.7) 


(interchanging i and j in the second sum). Interchanging the order of sub- 
scripts gives 


m m 
V 


j= 


'2 ifel ifel ' ^ 3' -X ' X , 3' 


where 


m m 



x=i j= 


11 1 11 

Tl--’.* T1 


(3.8) 


= L . ^ (M -M. ) 
' , 1 1 J 


= 6^ - 


i, j — • • • #ni 



- 10 - 


where the computation of 


1 D 


(3.9) 


involves a forward substitution, that is, is obtained from 




= m(3) A JU) 6 (3) 
P p ^ pq q 


p— 2,«..,n . (3.10) 


We thus have, from (3.8) that 



(3.11) 


where 

(ij ) (ii) (ij ) 
ri = 6 - 6 

P P P 

r (i j ) T 

and the j are calculated from (3.10). 

The use of the above formulae should probably not be compared with the 


approach used in LARSYS itself, but with the improvements proposed by 
6. Austin^^^ which take full advantage of the symmetry of the and of 

the symmetric structure of the summands in (3.7). It can be shown that the 
amount of v;ork involved in calculating in (3.11) is comparable with that 
involved with the corresponding terms in Ref (8). However, the amount of work 
involved in evaluating (3.4) is actually considerably less than the method 


-n- 


proposed in Ref (8) on account of the asymptotic linear dependence on m, 
as opposed to the quadratic dependence of Ref (8). It should nevertheless 
be pointed out that, from (3.4) 




'• 2 
= tr (KK) - nm 

(3.12) 

where 


m 




> 

II 

II 

1 

!-• 

(3.13) 


-1 




Thus, if the have been precomputed, the amount of work involved in 
evaluating may become negligible compared with the evaluation of 
by using (3.12) and the fact that, for symmetric matrices A, B: 


‘La 


a . .b . . 
ID Di 


n 


r" 



i-1 

\ . 

3^ 


a . .b . . 
ID Di 


+ 


n 




a . .b . . 

11 11 


However, this approach does not obviate the overall savings in feature selection 
of using the Cholesky decomposition instead of computing matrix inverses. 
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4. Classification : 

Classification involves the calculation of the maxiiinum likelihood 
functi ons 


fj(x) = cCj exp[- %(x-M)*Kj ^(x-M)] 

. j=l , . . . , 


(4.1) 


where x is the observation vector, = 1/(2 tt)®^^ , 


and 


-l.h 


CL. = (det K. ) 
3 3 


= in 


p=i 


(4.2) 


Actually, since exp(x) is a monotonic increasing function, only 

log f,. (x) needs to be computed in determining the maximum of f. (x) over 
J J 

all ‘ m classes. 

However, (4.1) is again simplified by noting that 

(x-M. ) * K. ^(x-M.) = y.* D. ^ y. 

3 3 3 -^3 3 ^3 


where 


y . = L . ^ (x-M.) 
•^3 3 3 


is calculated in a manner analogous to (3.9). 



-13- 


5. Results : 

The above techniques have been tested by appropriately modifying the 
OS version of LARSYS^ ^ supplied by NASA-JSC. These modifications are listed 
in Appendices A and B. In actuality the modification to the divergence cal- 
culations in feature selection use the Cholesky decomposition as opposed to 
the modified Cholesky decomposition as discussed in Section 3 - further 
savings of time, obtained by not having to calculate square roots, could be 
realized by using the modified Cholesky decompsosition. 

The modifications were written in single-precision FORTRAN and compared 
with the original single-precision versions in LARSYS. In the case of classi- 
fication, the results were also compared with a single-precision version of 
the corresponding calculations in LARSYS written in assembly code. 

The precision of these timing results is very open to question due to 
the difficulty of obtaining accurate and reproducible timing information under 
the OS Operating System of the IBM 370/155. Timings are heavily dependent on 
general system activity; furthermore the considerable subroutine overhead in- 
herent to the computation tends to mask much of the potential arithmetic 
economies of efficiency. 

The results are summarized in Figures 1, 2 and 3 on test data supplied 
by Purdue University with LARSYS. Figure 1, depicts the ratio of the time 
taken by the, original LARSYS version (DIVERG) to that taken by the proposed 
algorithm (CHOLESKY) in a divergence calculation for feature selection 
using six channels; this ratio is plotted for a varying number of target 
classes. It can be seen that CHOLESKY is approximately twice as fast as DIVERG. 
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Theoretical analysis shows that this ratio should greater than three 
for all values of m, and asymptotically should approatfn four for large values 
of m . This discrepancy underscores the high degree of tmprecision associated 
with the timing results. 

In Figure 2, the same ratio is plotted for a fixed number of classes (11) 
and where a varying number of features is selected from twelve channels. Ex- 
cept for a very small number of features, v;here the order of the is so 
small that the time of calculation is dominated by computaitional overhead, it 
can again be seen that CHOLESKY is between two and three times faster than 
DIVERG. Again, theoretically, this ratio should be betv^en three and four for 
all values of the number of features. 

In Figure 3, the time taken for classification using the three methods 
is compared for a number of points varying from 50,000 to 100,000. The 
Cholesky method is significantly faster (about 30%). It should be pointed out 
that, as has been noted el.sewhere^^^ , equivalent savings can be obtained by 
using a variant of the LARSYS calculations which does not employ the modified 
Cholesky decomposition; however, this variant does not have the accuracy 
potential of the Cholesky approach^^). 

6. Improvements in Accuracy : 

The modifications described were executed in single-precision so as to 
provide a basis for comparison with the LARSYS calculations. Without further 
refinement, it should not be surprising that the accuracy will be corres- 
pondingly limited, since^^^ accuracy in such computations is essentially a 
function of three principal components: 
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the method employed 
. the arithmetic significance 

. the conditioning of the various matrices 

For ill-conditioned systems (in the applications unfe consideration, 
these may arise, for example, from working with highly- correlated channels), 
more precise methods have to be -employed and/or the aritteietic significance 
increased. Directions which need to be examined with higlhier accuracy ob- 
jectives in mind include, not only that of using higher significance arith- 
metic, in sensitive portions of the computation, but also those of employing 
iterative refinement, scaling or equilibration. These will;, however, be 
studied in a later report. 
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Figure 1 

Timing comparison as a function of number of classes used for the 
divergence calculation. 





Figure 2 

Timing comparison as a function of number of channels selected for the divergence calculation. 





APPENDIX A 

Flowcharts and Listings of Subroutines 
Used in Performing Feature Selection 
Using the Divergence Criterion Employing 
the Gholesky decomposition 


DIVGAL 


Subroutine DIVGAL, which performs the divergence calculation, and 
stores the channels corresponding to the maximum separation, is called by 
the main program, which reads the aircraft data tape and calculates the mean 
vectors and covariance matrices of the various classes from specified 
training fields. 

On the initial call to DIVGAL, the mean vectors and covariance 
matrices of each class corresponding to channels 1 to NF, where NF is the 
number of features desired, are placed into vectors in DIVGAL, The com- 
putation is done in vector rather than matrix form, since the computer will 
handle calculations with single subscripts faster than those involving 
multiple subscripts. However, the expressions which follow are given in 
matrix form for the sake of clarity. 

The search for the set of channels which exhibits the maximum sepa- 
ration is done by the exhaustive technique in which every possible combina- 
tion of NF channels is examined. Because of the method used to step 
through all the combinations, often the first (NF-1) channels selected remain 
unchanged, with only the highest numbered channel changing between diver- 
gence calculations. Thus, in subsequent calls, only those elements of the 
mean vector and covariance matrices corresponding to the changing 
channels are redefined, with the other elements being retained in storage. 

DIVGAL then performs the Gholesky decomposition of the covariance 
matrix for each class, and simultaneously accumulates the sum of the 
matrices. The Gholesky decomposition of the sum of the covariance 
matrices is then calculated. 
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The total interclass separation is computed from 


TOTSEP = D1 + D2 

where D 1 and D2 are function subprograms. D 1 is defined as 

Dl=I||L'ljL iP 

where the sum over i is over all the classes, L ^ is the inverse of the 

j 

J fVi 

Choleskj'’ decomposition of the NF X NF covariance matrix of the j 
class whose elements correspond to the set of channels under consideration, 
L is the Cholesky decomposition of the NF X NF sum of covariance 
matrices likewise corresponding to the set of channels under consideration 
and 1 1 A I p is defined as 

l|A|P=Z 

ij iJ 

0 

D 1 performs H A jp and the sum over classes after calling subroutine CK 
to perform LT^^L by back substitution 

D2 is defined as 

D2 = ^ i j 

1,3 

where r].. = , 6.. = L IVL, and M. is the mean vector formed of 

ij 11 ij ij 1 J ’ J 

the selected channels of the class. D2 performs the multiplication and 
sum over classes after calling subroutine DJM to calculate by back 

substitution. 
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SUBROUTINE DIVCAL 
















2 



•4 

Jconvert matrix of covariance 
.matrices into vector 

i 


iperform Gholesky decomposition 
:of covariance matrix of 
■ KLS^“ class 













■6 

: IMSL subroutine LUDEGP returns 
j diagonals in reciprocal form- thus 
j_ recipro^cal of diagonals taken* 


I 7 ( 

I / ( 

/■! Accumulate sum of covariance 
: matrices 1 

I — j 


* IMSL , 6200 Hillcroft , Suite 510 , Houston , Texas 77036 











: perform Gholesky , 
; aecomposition of i 
; sum of covariance i 
1 matrices ! 
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D, (L, SUML, NF, NUMCLS, LENGTH) 



, C3-Icul9.t0 Lj JCL-/S ^ 
1 by back substitution 


j calculate I jL 1 ^tid [ 

[accumulate sum ove_r_classes| 


"'calculate D1 ; 
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CK(L,SUML,C,NF) 



.-9 





































D2 (L, SM, NF, NUMCLS, LENGTH). 



DEE2=0 



1 ~ 1 ‘ 
I calculate Lj*m. ; 

Iby back substitution 


i -1 

: Calculate L ^ m. 

: * 

I by back substitution ; 
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NOV 71 ) 


OS/360 FORTRAN H 


;OMPILER OPTIONS - NAME= M A I N , O? T= 0 2 . L 1 NEC NT = 6 0 t SI ZE = 0 00 OK , 

SOURCE . EBCDI C. NOL I ST ,NOOECK,LOAO, MAP, NOE DI T. I D, XREF 
INTEGER’f‘2 BLOC K( 4 ) , C SE L ( 1 2 ) / 1 2 *1 BD AT A ( 1 500 ) 

IN T EGER* 4 ID{ 200 ) , F R ( 5 , I 2 ) . ERROR , C SET ( 3, 1 2 ) , NN H 6, 30 , 1 0) ,NN2( 30) ♦ 
ICH SEL( 12 ) , BSTCHI 12 ) 

REAL *4 RDA TA( 3000 ) , SDATAI 1 2» 10 00) , DATTOT ( 1 2) .TU 2 , 1 2 ) .MEAN! 12,30 } , 
1 K( 1 2 , 1 2 , 30 ) ,S1GM A( 1 2 , 30 ) , RHO ( 1 2, 1 2, 30 ), TOT SMP{ 1 2 ), MA XSEP 

* REALMS CLSl 30) 

EQUIVALENCE ( I D{ 51 ) , FR { 1 , 1 ) ) 

COMMON C LS , MEAN, S I GM A, K , RHO , FR , N N 1 , NN 2 , T OT S MP 

ID ( 1) =0 

ID(2)=0 

NCO=12 

CALL GADRUNI 66000600 ,3 , ID , ERROR) 

IF { ERROR ,GT .0 ) GO TO 35 

DO 10 1=1 ,3 

DO 10 J= 1 , 1 2 

CSET ( I , J ) = FR( 1+2 , J ) 

10 CONTINUE 

READ NUMBER OF SETS OF TRAINING FIELDS ' 

C . 

READ { 5, 15) NUMSET 
DO 150 NUMRUN= 1 , NJMS ET 
C 

C‘ . READ TOTAL NUMBER OF CLASSES 
C 

READ (5,15) NUMCLS 
1 5 FORMAT! I 5, 1 X, A 8) 

DO 80 KLS= 1 , NJMCLS 
NSMP=0 
C , 

C READ NUMBER OF TRAINING FIELDS (NN2(KLS)) IN CLASS CLSIKLS) 

C . ■ 

READ (5,15) NN2( KLS) ,CLS(KLS) 

C \ ■ ■ 

C INITIALIZE ACCUMULATORS FOR UNNORMALIZED MEAN VECTOR 

C .. AND COVARIANCE MATRIX .. . ... 

C-- 

DO 20 1=1 , 1 2 
DATTOT{ I )= 0. 

DO 20 J= 1 , 1 2 
T( I , J)=0. 

20 CONTINUE . _ . 

NN2 KLS -NN2 ( KLS ) 

DO 65 I1=1.NN2KLS 
C 

C READ BOUNDARIES OF Il(TH) TRAINING FIELD OF THE KL S( TH ) CLASS 

C NNl ( 1 ) = FIRST LINE NNI(4) = COLUMN INCREMENT 

C NNK2) = FIRST COLUMN NNK5) = LAST LINE 

C NNl (3) = LAST COLUMN NNK6) = LINE INCREMENT 

C 

* READ ( 5, 25 ){ NNl ( I , KL S, 1 1 ) . I = 1 , 6) 

25 FORMAT (6 15 ) 

C 

. C COMPUTE number OF SAMPLES PER CHANNEL PER. LINE 

C 


N1 =NN1 ( 1 ,KLS,I 1 ) 
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oo nno onn onon 


N2=NN1 (2 «KLS,11 ) 

N3=NN 1 ( 3 ,KLS, I 1 ) 

N4 = NNJ (4 ,KLS, I 1 ) 

N5 = NN1 ( 5 «KLS ,I 1 ) 

N6 = NN 1 ( 6 ,KLS, I 1) 

NSD = ( ( N3-N2 )/N4 ) +1 
NRQ=NSD+6 
C 

C PREPARE BLCCK(I) FOR G ADL I N 

C 

00 30 I=2t 4 

BL0CK( I ) =NN1 ( I ,KLS , I 1 ) 

30 CONTINUF 

READ IN DATA FROM I2(TH) LINE OF IHTH) TRAINING FIELD 
OF ■ -THE ’KUST-Th.) .CLASS 

DO 60 I2=N1 , N5 ,N6 

BLOCK! 1 ) = I 2 

CALL G ADL I N( BLOCK, CSEL , CSET, I D . 3 , NCD , NR 0 , BD A T A , RD A TA , ROLL .ERROR) 
IF ( ERROR .GT .0 > GO TO 00 

SAVE THE DATA 

DO 45 1=1,12 
KK2 = ( I-l )»''NRQ 
DO 45 J= 1 , NSD 
SDATA( I , J) =RDATA ( J+KK2 ) 

4 5 CONT INUE 

ACCUMULATE MEAN VECTOR 

00 50 1= 1 , 1 2 ■ 

DO SO J=1 , NSD 

DATTOT ( I )= DATTOT ( I ) +SDATA( I , J ) . 

50 CONTINUE 

ACCUMULATE COVARIANCE MATRIX 

DO 55 1= 1 , 1 2 
DO 55 J= 1 , 12 
DO 55 L=I,NSD 

T( I , J ) =T( I , J) + SDATAI I , L) *SDATA ( J ,L) 

5 5 CONT INUE 
NS MP=NSMP+ NSD 

60 CONTINUE 
55 continue 

TOTSMPI KLS)=NSMP 
C 

C . CIV IDE ACCUMULATORS BV TOTAL NUMBER OF SAMP_ES TD 03TAIN 
C MEAN VECTOR AND COVARIANCE MATRIX OF KLS(TH) CLASS 

C ' 

DO 70 1= I , I 2 

MEAN! I ,KLS )=DA TTOTd ) /TOTSMP (KLS. ) . . 

7 0 CONTINUE 

DO 7 1 1 = I , 1 2 , 

DO 71 J= I, I 2 

K( 1 , J , kLS ) =T{ I , J ) /TOTSMP! KLS) -ME AN( I , KLS) *MEAN( J ,KLS) 
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K( J» I,KLS)=K( I , J ,KLS) 

71 CONTINUE 

DO 75 I = 1 » 1 2 

S IGMA( I.KLS) = SQRT(K{ I ♦ I ,KLS) ) 

DO 75 J= 1 . I 

RHO( I . Jf KLS) =K( I , J .KLS) /SORT (K( I ♦ I .KLS ) *K( J, J, KLS) ) 

RHO ( J, I . KLS ) =RHO { I , J , KLS ) 

7 5 CONTINUE 
80 CONTINUE 

CALL ST ATPT (NUMCLS ) 

GO TO 100 

85 CALL E RSDRN{ ERROR) 

GO TO 955 

90 CALL ERGDLN(ERROR) 

GO TO 955 
100 -CONTINUE 
C 

C READ NUMBER OF FEATURES TO BE USED IN SELECTION 

C 

10 5 READ ( 5*25) NF 

IF (NF.EQ.O )GQ TO 95 
WRITE (6,1000) NF 

1000 FORMAT {• IFEATURE SELECTION FOR*. 15,* CHANNELS*) 

c 

C initialize CHANNEL SELECTORS 

C- 

DO I 1 0 1 = 1 , NF 
CHSEL(I)=I 

110 continue 

LC=NF 
KOUNTR = 0 
11 = 1 
CD=NCD 
FE.AT=NF 
MAXSEP=0. 

NCC = 1 

LENGTH = NF« (NF+l ) /2 
C 

CALCULATE NUMBER OF DIFFERENT COMBINATIONS TO BE TESTED 
NUMCHK = GAM MA(CO+l.)/(GAMMA( = EAT-H.)>!'GAMMA(CD-FEAT+l.)) + .l 
PERFORM DIVERGENCE AND SEPARABILITY CALCULATIONS 
CALL RT I ME ( JBEGIN ) 

115 CALL DIVCA L{ K, MEAN, NUMCLS .NF .LENGT H.CHSEL. BSTCH, M AXSEP, NCC) 
KOUNTR=KOUNTR+ 1 

\ 

^CC = NF 

CHECK WHETHER ALL POSSIBLE COMBINATIONS HAVE BEEN EXAMINED 

IF { KOUNTR. EQ.NUMCHK ) GO TO 135 
CHSEL(LC) = CHSEL( LO + l 
1F(CHSEL(LC) .LE.NCD)GO TO 115 
120 NN=LC-II 

CHSEL( NN ) = CHSEL( NN) 4- 1 
KCC=NN 

IF{ CHSEL( NN) .G T. ( NCD-I I ) ) GO TO 130 
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NR = NN+1 

DO 125 JJ=NR,NF 
ChSEL( JJ )=CHSEL( JJ-1 )+ 1 
125 CONTINUE 

11 = 1 

GO TO 1 1 5 
13011=11+1 

GO TO 120 

135 CALL RTIME(JEND) 

WR I TE( 6, 95 95 ) ( CL SI KL S) ,KLS = 1 ,NUMCLS) 

9595 FORMAT {* ICLASSES: */5X, 10(A8« 3X)/5X, 10(A8,3X)/5X. 10(A8»3X)) 

W.R I TE (6 ,9 000) NF 

9000 FORMAT {* OFEA TORE SELECTION FOR*. 15, • CHANNELS') 

WRITE 15, lAO) NF, (BSTCHI r ) , 1= 1 ,NF) 

140 FORMAT! 'OTHE BEST* ,15,' CHANNEL(S) ARE', 1215) 

-EL AP-T-M = ( dB BGtl'N-dENDO /I 00, 

WRITE 16 ,200 ) NUMCLS , NF, EL APTM 

200 FORMAT!//' ELAPSED TIME FOR FEATURE SELECTION FOR', 15, • CLASSES,' 
1,15,' CHANNELIS), I S ‘ , F 8 , 2, • SECOND S ' ) 

WR I TE ! 5 ,1 41 ) 

14 1 FORMAT! IH 1 ) 

GO T O 1 0 5 

95 CONTINUE 
150 CONTINUE 
955 ST OP 

End , ' 
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NOV 71 ) 


OS/360 FORTRAN H 


compiler options - NAME= M AI N , OP T = 02 , L INEC NT =60 ♦ SI ZE =0 00 OK. , 

SOURCE .EBCDI C. NOLI STf NODECK, LOAD, map .NOEDI T t lO.XREF 
SUBROUTINE S TA TP T( NUMC LS) 

INTEGER*4 FP { 5 , 1 2 ) , N N I (6 , 30 , 10 ) , NN 2{ 30 ) 

RE AL*4 MEAN( 12 ,30 ) ,S IGMA ( 1 2. 30 ) . K( 12 , 12, 30 ) ,RHO( 12, 12, 30 ) , 

ITOTSMP (12) 

REALMS CLS(30) 

COMMON C LS , MEA N, SI GM A, K , RHO, FR ,NN1 , NN2 , T OTSMP 
' DO 10 KLS=1,NUMCL5 

WRITE (6,100) KLS, CLS (KLS ) 

100 FORMAT! • IMEAN VECTOR, COVARIANCE MATRIX, AND CORRELATION MATRIX FO 
IR CLASS NUMBER*, I 3,* (*,A8,*)*/) 

NTF=NN2(KLS) 

DO 35 11=1 ,NTF 

i WRITE (6,05) I 1 , NNl ( 1, KLS, i 1 ) ,NN 1( 5, KLS , I 1 ) , NNl ( 6, KL S, I I ) ,NN 1( 2 , KL 

'1'S,'I'10', NMH 3 ,KLS-, II-)-, NNl (4, KLS, M ) 


95 FORMAT! 5X, *TRA INING FIELD NUMB ER * , I 4 , • : LINES’, 16,* TO*, 16,* (EVER 

1Y*,T3,* LINE(S)), SA MPLES * , I 6, * TO', 16,’ (EVERY*, 13,* SAMPLE(S))*) 

35 CONTINUE 

WRITE ( 6 , 9 6 ) TOTSMP! KL S) 

96 FORMAT ( I OX , * TOT AL NUMBER OF DATA POINTS = •,F7.0//> 

WRITE (6,110) { FR( 1 , I ) , I =1 ,1 2 ) , (FR (2 , I > , 1 = 1 , 12 ) 

110 FORMAT ( 1 X, *SPECTRAL *,12(F5.2,*- * ) /3 X , * BA NO • , 4X , 1 2 { F 5 , 2 , 4X) / ) 

WRITE (6,120) ( MEAN! I, KLS) , 1 = 1 ,1 2) , ( S IGMA! r,KLS) , 1 = 1 , 12) 

12 0 FORMAT { * MEAN* , 5X, 1 2 ( F 7.2 , 2X) //* ST DEV*,3X,12(F7,2,2X)) 

WR ITE (6 , 1 25 ) 

125 format!///* COVARIANCE MATRIX*/) 

WRITE (6,110) ( - R( 1, I ) ,I = 1 , 1 2) ,( FP ( 2. , I ) , 1 = 1 , 12 ) 

DO 20 1 = I , 12 

WRITE (6,130) - R ( 1 , 1 ) , (K( I , J ,KLS ) , J=] , I ) 

130 FORMAT ( 1XF5. 2, * - • , 2X, 1 2(F 7.2,2X) ) 

WRI TE (6 , 1 40 ) FR (2 , I J 
140 FORMAT ( 2X, F5. 2) 

20 CONT INUE 

WRITE (6,126) 

126 FORMAT! ///• CORRELATION MATRIX*/) 

WRITE (6,110) { FR( 1 , I ) , 1= 1 , 1 2) , ( FR ( 2, I ) , 1= 1, 12) 

00251=1,12 

WRITE (6,130) FR ( 1 , I ) , (RHO ( I , J ,KLS ) , J=I , 1 ) 

WR I TE (5 , 1 40 ) FR ( 2 , I ) 

25 CONTINUE 
10 CONTINUE 
RETURN 

end 
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NOV 71.) 


OS/360 FORTRAN H 


COMPILER OPTIONS - NA ME= M A I N « OPT-0 2 . L I NEC NT =6 0 . S I ZE=0 0 0 0 K , 

SOURCE .E BCD I C , NOLI ST.NODECK, LOAD, MAP, NOE OIT,ID,XREF 
: SUBROUTINE ER GD R N ( E RRO R ) 

! INTEGER#4 ERROR 

WR ITE ( 6, 1 0 ) ERROR 

i GO TO ( 1 01 , 1 02 , 1 03, 1 04 , 105, 106 , 1 07, I 08) , ERROR 

, > 10 1 WR I TE ( 6 , 1 ) 

GO TO 109 , 

i - 102 WRI TE (6 ,2 ) 

' GO TO 109 

■ I 103 WRITE (6,3) 

GO TO 109 

: 10 4 WR I TE ( 5, 4 ) 

; GO TO 10 9 

'■ 1 0 5 WR I TE { 5 , 5 ) 

I GO TO 10 9 

. 10 6 WR ITE (6 ,6 ) 

GO TO 109 
10 7 WR ITE (6,7) 

I GO TO 109 

■I 1 0 a WR I TE ( 6 , a ) 

10 9 WR ITE (6,9 ) 

1 FORMAT (5 X, * EXPECTED ID RECORD READ AS END OF FILE RECORD*) 

1 2 FORM AT { 5X, ■ EXPEC TED ID RECORD HAD WRONG READ COUNT*) 

3 FORM AT (5X, • EXPECTED ID RECORD READ WITH A TAPE P AR I TY CHECK « ) 

'i 4 FORMAT! 5X, • E XPEC TED ID RECORD READ WITH A HARDWARE PARITY CHECK*) 

1 5 FORMAT! 5X, ‘EXPECTED ID RECORD HAD WRONG READ COUNT AND PARITY CHEC 

TK* ) 

r 6 FORMAT! 5X ,* TAPE UNIT NOT ASSIGNED* > 

i 7 FORMAT ( 5X, *EXPEC TED ID RECORD READ FROM CORRECT TAPE AND FILE, BUT 

1 CONTAINED WRONG RUN NUMBER*) 

> 8 FORMAT! 5X, • RUN WAS NOT FOUND IN ADCRT*) 

9 FORMAT (/ 10 X, 'EXECUT ION TERMINATED*) 

'■ 10 FORMAT (• lERROR NUMBER*, 13.* IN SUBROUTINE GAORUN •/ ) 

RETURN 

.ENTRY ERGDLN (ERROR ) 

; WRITE (6,28) ERROR 

j GO TO (1 1 1,1 12, 1 13,1 14,1 15,1 16,1 17,1 18, 1 19,120,121 ,122,123,124,125 

1 ,1 26 , 1 27 ) , ERROR 
I 1 1 WR I TE ( 6 ,1 1 ) 

GO TO 128 
1 1 2 WR I TE (5 , 1 2 ) 

■ GO TO 128 

•11 3 WR ITE (6, 1 3) 

GO TO 1 28 

1 1 4 WR I TE ( 6, 1 4) 

' , GO TO 128 

1 I 5 WR I TE ( 6 , 1 5) 

‘ GO TO 128 

' 1 1 6 WR ITE (6,16). _ 

GO TO 128 V 

1 17 WR ITE (6,17) 

' GO TO 128 

1 1 8 WR I TE ( 6 , 1 8) 

GO TO 128 • ■ ' 

119 WRITE (6,19) 

' GO TO L2 8 

120 WR ITE (6,20 ) 
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GO TO 128 
12 1 WR ITE (6,2 1) 

GO TO 128 
122 WRITE (6,22) 

GO TO 128 
12 3 WR I TE ( 6, 2 3) 

GO TO 128 

124 WRI TE (6 ,24 ) 

GO TO 128 

125 WRITE (5,25) 

GO TO 128 

126 WR I TE ( 6, 2 6 ) 

GO TO 128 

1 2 7 WR I TE ( 6 , 2 7 ) 

12 8 WR ITE (6,9) 

1 1 FORMAT (5X , ’DAT A LINE REQUESTED DOES "NOT EXIST ON TAPE') 

12 FORMAT! 5X, ’BYTE COUNT IN REQUESTED DATA RECORD IS INCORRECT*') 

13 FORMAT ( 5X ,* PAR IT Y CHECK OCCURED IN READING REQUESTED DATA RECORD') 

14 format (5X HARDWARE ERROR CCCURED IN READING REQUESTED DATA RECORD 
1 » ) 

15 FORM AT (5X , • IN READING REQUESTED DATA RECORD COMBINATION ERROR OF • 
1 /5X’ I NCORRECT BYTE COUNT AND PARITY CHECK DR • / 5 X * I NCORR ECT BYTE CO 
2UNT AND HARDWARE ERROR') 

16 FORMAT (5X ,' T APE UNIT NOT ASSIGNED') 

17 FORMAT! 5X LINE WAS DEFINED LESS THAN OR EQUAL TO ZERO') 

18 FORMAT (5X, * I SAM WAS DEFINED LESS THAN OR EQUAL TO ZERO') 

19 FORMAT (5X , 'S INT WAS DEFINED LESS THAN OR EQUAL TO ZERO') 

20 F0RMAT(5X, ' I SAM IS GREATER THAN LSAM') 

2 1 FORM AT (5X , ' CSEL FLAGS ARE LESS THAN 0 OR GREATER THAN 7') 

22 FORMAT (5X NO CHANNELS SELECTED OR SELECTED CHANNELS NOT IN RUN') 

23 FORMAT ( 5X, 'NCR IS GREATER THAN NCD • ) 

2 4 FORMAT (5 X, 'NSO IS LESS THAN NSR + 6*) 

25 FORMAT ( 5 X, 'DATA IN REQUESTED LINE. DOES NOT EXIST') 

26 FORMAT (5X , 'DATA CANNOT BE CALIBRATED AS REQUESTED' > 

27 FORMAT (5 X, • REQUESTED LINE CANNOT BE LOCATED ON TAPE*) 

28 FORMAT (• IE RRCR NUMBER', 13,' IN SUBROUTINE GAOLIN'/) 

RETURN 

END 
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NQV 7 1 ) 


OS/360 FORTRAN H 


COMPILER OPTIONS - NAME= M A I N * OPT = 0 2 . L I N ECNT = 6 0 . S 1 2 E= 0 0 0 OK , 

SOURCE »EBCDI C.NOLI ST , NOD ECK, LO AD , M AP , NOE D I T , ID,XREF 
SUBROUTINE D I V C A L ( K , ME AN , N UM CL S , NF »L ENGT H , C H SE L , BS TC H , M AX SE P ,N CC ) 
RE AL*4 K ( 1 2 , 12 , 1 ) , ME AN ( 1 2 ♦ 1 ) , L < 2 34 0 ) . SM { 36 0 ) , S K ( 23 40 ) . SU MK ( 7 8) 
REAL*4 MAXSEP, SUML( 78) 

INTEGER*4 C HSEL ( 1 2 ) , 8S TC H( 1 2 ) 

DO 5 1=1 .LENGTH 
5 SUMK( I )= 0. 

NSK=-LENGTH 
! NSM=-NF 

DO 20 KLS=1,NUMCLS 

NS K = NSK+ LE NGTH 

NSM = NSM+NF ‘ 

II=NSM+NCC-1 
I I I =NS K+ NCC« (NCC-1 ) /2 
DO 10 T=NCC-.’NF 
I 1= I I + l 
C 

C CONVERT- MA TRI X OF MEANS INTO VECTOR 

C 

SM( I I )=MeAN (CHSEL ( I ) .KLS ) 

; . DO 10 J= 1 , I 

I I I = I I I 1 
C 

CONVERT MATRIX OF COVARIANCE MATRICES INTO VECTOR 

10 SK( I I I )=K{ CHSEL ( I ) . CHSELI J ). KLS ) 

NSPl=NSK-t-l 
IER=0 

PERFORM CHOLESKY DECOMPOSITION OF MATRIX 

CALL LUDEC P ( SK (NSPl ) .L (N SP 1 ) , NF , DD 1 ,DD2, lER ) 

• I F ( I ER .NE . 0 )RE TURN 

.DIAGONALS RETURNED BY LUOECP ARE RECIPROCAL 

ND=NSK .. 

DO 1 5 J= 1 . NF 
ND=ND+J 

15 L( ND)=1 ./L (ND) 

. ACCUMULATE SUM OF K MATRICES 


NL=NSK 

DO 17 1=1, LENGTH 

NL = NL-H 

17 SUMK{ I ) =SUMK( I )-! SK (NL) 

20 CONTINUE 

I 

COMPUTE CHOLESKY DECOMPOSITION OF SUM OF COVARIANCE MATRICES 

CALL LUOECP (SUMK ,SUML ,NF, ODD 1. ODD2, lER ) 

IF(IER.NE.O) RETURN 

DIAGONALS RETURNED BY LUDECP AR E R EC I PROC AL 
NDIAG=0 
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00 25 1 = 1 . NF 
ND I AG=N0 1AG+ I 

25 SU KLl NDI AG )=1 . /SUML (NOI AG) 

COMPUTE TOTAL IMTERCLASS SE^A^ATION AND COMPARE WITH MAXIMUM 
C SEPARATION OBSERVED SO FAR 

C 

TOTSEP=Di ( L , SUML . NF, NUMCLS .L ENGT H) +D 2< Lt SM, NF, NUMCL S, LENGTH ) 

IF ( TOTSEP. LE.MAXSEP) RETURN 
MAXSEP=T OT SEP 

WRITE (6,151) MAXSEP, ( CHSEL{ I ) , 1 = 1 ,NF ) 

151 FORMAT! • 0**=f==f‘TOTAL SEPARATION = ‘.IPEIA.S,* CH ANNELS : • , 1 21 5 

DO 50 1 = 1 , KF 

50 BSTCH ( 1 ) =CHSEL ( I ) 

RETURN 

END 


I 


A-23 



o o 


NOV 7 1 ) 


OS/360 FORTRAN H 


compiler options - NAME= M A I N , OPT= 0 2 » L I N ECNT = 6 0 , S I Z E = 0 00 OK , 

SOURCE ,E8CDI C, NOLI ST , NODECK, LOAD, MAP, NOEDIT, I D, XR EF 
REAL FUNCTION O 1 ( L , 5 UML , NF , N UM CL S , L ENG TH ) 

C 

C CALCULATION OF SUM OVER J OF jjLlJ) INVERSE # SUM Lli«*2 

C BY EACK SUBSTITUTION 

C 

PEAL=«=4 L( 1 ) , SUML ( 1 ) , C(78) , SM ( 1 ) , DI I( 1 2 ) ,D( 1 2 ) 

' SUMC=0. 

NS K = 1 -LENGTH 

■ ^ DO 20 KLS=1,NUMCLS 

NS K-NSKFLENGTH 
CALL CK( L( NSK ) , SUML , C, NF ) 

DO 10 11=1 .LENGTH 

Cl l-Cl 11 ) 

10 SUMC. = SUMC+'CIT*CT I 
20 CONTINUE 

Dl =SUMC-NF’!=NUMCLS’I=NJMCLS 
RETURN 

entry D2 ( L ,SM, NF , NUMCL S, length) 

C 

CALCULATION OF SUM OVER I AND J OF ET A { I , J ) ( * ) =i=ET A ( I , J > 
WHERE ET A( I, J ) = OElTA ( I , I >-DELTA( I , J ) , AND 
C DELTA( I , J) =L( I) INVERSE^M ( J) , 

C BY BACK SUBSTITUTION 

C ... . 

DEE2=0, , 

LN SM= 1-LENGTH 
LNSV=1-NF 

DO 100 LKLS-=1 , NUMCLS 
LNSM=LNSM+ LENGTH 
LNSV = LNSV + NF , 

CALL D JL( L { LNSM) ,SM( LNSV) , DI I , NF ) 

JNSV=1-NF 

■ DO 100 JKLS=1 , NJ MCLS 
. JNSV= JNSV+ NF 

IFILKLS.EQ .JKLS ) GO TO 100 
CALL DJL(L (LNSM) ,SM( JNSV ), D,NF ) 

DO 1 20 I -1 ,NF 
ET A=D I I( I )-D (I ) 

120 DEE2=OEE2 + ET A*ET A 
lOOCONTINUE 
D2=DEE2 

RETURN ' . 

end 



O U 


NOV 71 ) 


OS/360 FORTRAN H 


COMPILER OPTIONS - NAME= M A IN , T= 0 2 . L IN EC NT =6 0 » S I ZE =0 00 OK , 

SOURCE, EBCOI C, NOLI ST , NOD EC K, LOAD , MAP , NOEO I T , ID, XREF 
SUBROUTINE CK( L , SUML ,C , NF ) 

C 

C SUBROUTINE TO CALCULATE L( J ) INVERSE * SUM L BY BACK SUBSTITUTION 

C 

RE AL*A L ( I ) , SUML { I ) , C( 1 > ,SM( 1 ) ,D { 1 ) 

c{ 1 ) -suML( 1 ) /L n ) 

' IF { NF .EQ . 1 ) RE TURN 

NBEGIN = 1 

/ 00 1 00 NRW=2 ,NF 

NBEG IN = NBE G I N + NR W- 1 
NEND = NB£GI N+NRiV- 1 
NE M1=NEND-1 
NSTCR-0 

DO -50 T-NBEGIN ,'NEMl 
S=0. . 

ICR^I-NBEGIN 
NSI CR=N5 ICR+ICR 
NC=NS1CR+1 
DO 20 NL=I,NEM1 
NC =NC+ ICR 
ICR=ICR+ 1 

20 S=S+L{NL )*C(NC) 

50 C( I ) =( SUML ( I )-S) /L (NEND) 

10 0 C( NEND )= SUMLINEND ) /L ( NEND) 

RETURN 

ENTRY 0 JL( L ,SM ,D ,NF ) 

SUBROUTINE TO CALCULATE L(I) INVERSE * M(J) BY BACK SUBSTITUTION 

D ( 1 ) = 5M( 1) /L ( 1 ) 

IF{ NF .EO.l IRETURN 
NBEGIN=1 

. DO 200 NRW = 2,NF 
■N0EG 1 M=NBEG I N+NRW-1 
NEND-NBEGI N+NRW-1 
NEM 1 =NEND- 1 , . _ 

S = 0. 

1 = 0 

DO 210 NL= NBEG IN , NEM I 

1 = 1 + 1 

21 0 S=S+L(NL >^D( I ) 

200 D(NRl¥ ) = ( SMINRW )-S)/L (NEND) 

return 

EN D 


I 

! 


A- 25 


/^^o^PoZT 



Subroutine MCHLSK 


This subroutine computes the modified Cholesky decomposition of the 
covariance matrix as described earlier. The decomposition overwrites the 
covariance matrix. The determinant of the covariance matrix is also cal- 
culated. The covariance matrix is stored in symmetric storage mode (i.e. - 
the upper triangular part is stored by columns). The diagonal elements of 

the diagonal matrix D are stored as - 1 in the diagonal positions of the 

2TJi i 

covariance matrix. The off-diagonal elements of the lower triangular matrix 
are stored in the corresponding positions in the covariance matrix (the dia- 
gonal elements of L are all equal to unity). To simplify the understanding 
of the flowchart, full matrix notation has been used for matrix elements. 
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START 



Det'=*— 1 


TF«-KK(J,J) "• 


J = 1 


i No 


Duni(I>— KK(I,I)*KK(J,I) 

TF - Dum (I)*KK (J,I) 
I+l 


KK(J,J)-5— TF 


DET*TF 


On input KK is the covariance matrix 
NV is the number of channels used. 
DET will be the determinant. 


• Loop to compute D (1,1) and store in 
TF, then transfer it to KK(J,J) 


j Save Dum (I) for 
• later use 
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Subroutine CLASS 


This subroutine classifies N0 data points using the modified Cholesky 
decomposition. The decomposition is generated by subroutiliire MCHLSK. The 
class number selected and the corresponding value of the dtaisity function are 
saved in arrays IR and VR, respectively. Arrays COR and iHM. are stored and 
indexed as vectors, but, for simplicity, we have used matriix notation in the 
flowchart. 
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START 


TFMAX<f— 

>_io60 

II <s-l 



RDATA contains the data for N0 vectors 
(stored by channels). NC is the number 
of classes. NV is the number of channels 
used. C0N contains the constants used to 
compute the density function. C0R contains 
the modified Cholesky decompositions 
(L(I,0) and - 1 ) . AVE contains the 


IBd^5 — II - N0 - 6 

— I 


mean vectors. 


IBD< — IBD + N0 + 6 
DATA (I)<^RDATA (IBD; 
I<s— I + 1 


-) Move elements of II^*^ vector into 


I Loop over all classes ! 

I I 


DATA(l) - AVE (JJ,1) 

DM(1)<5 — S 

TF<^C0N(JJ) + C0R(JJ,1,1) *S*S 


KD 2 


Loop for computing the KD 
Y; store it in S 


element of 


0ATA (KD) - AVE (JJ,KD) 
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Subroutines COVIN AND STATS 


These LARSAA routines were slightly modified to employ the modified 
Cholesky decompositions. COVIN now calls MCHLSK and does not compute K“^. 
STATS now employs MCHLSK to calculate the determinants. 
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APPENDIX C 


Flowcharts of Routines Necessary 
to Employ the Modified Cholesky 
Decomposition in Divergence 
Calculations 



Modified GHOLESKY Feature Selection Method. 


The accompanying flow charts detail the method of feature selection 
using the modified Cholesky decomposition, DIVCAL (Appendix A, pp. A-1 to A- 7) 
will be modified in that the call to LUDECP will be replaced with a call to 
a routine which will decompose K into LDL* instead of LL*. This may be 
accomplished by a subroutine such as MCIdLSK (see Appendix B), returning 
from the routine prior to the modification of the diagonal elements (i. e, , 
before the step preceding off -page connector D on page B-'3). Also, the steps 
in DIVCAL obtaining the reciprocal of the diagonal elements of L, returned 
by LUDECP, will be removed. The form of D1 and D2 will be changed, 
although the structure will remain quite similar. ; 

Although the summations are given in matrix form for reasons of 
clarity, the calculations should be performed using singly- subscripted 
variables. 

This modification has not yet been tested but the logic, as evident 
from the flowcharts, almost parallels that of the unmodified Cholesky method. 

For simplicity, and to save storage, the elements of and D are stored as 
the diagonal elements of Lj and L respectively. 


C-1 



D1(L,SUML,NF,NUMCLS, LENGTH) 





SUMC 
NSK = 

, 

= 0 

1- LENGITI 


2 . 

KLS = 1 


NSK-NSK+ LENGTH 


1 

CALL TKL(NSK), SUML, T, NF) 



I SU MC=SUMC+T(1) »T(1) *SUML(1)/L(NSK) 


W3— 




NBEGIN = 1 
NRW = 2 


' -1 ' ! 

j calculate L | 

:by back substitution ; 


NBEGIN = NBEGIN+NRW-l 
NEND = NBEGIN+NRW-l 
NDJ=: NSK+ NEND-]_ 


I = NBEGIN 


N1>1 
II = 1 



HSUMC = SUMC +T(I):icT(I)*SUML(ND)/L(NDJ)| - 


accumulate sum over i 
; elements and classes ofj 

i (tj ^ /dj ) d • ; 

' ^ pq^ q : i 



C-2 












TJ(L,SUML,T NF) 



[ Ta)=l 



C-4 






















D2 (L, SM, NF, NUMGLS, LENGTH) 


DEE2-0 


LNSM= 1- LENGTH 
LNSV = 1- NF 


[LKL S =_1 


rLNSM'= LNSA4+ LENGTH' 
LNSV = LNSV + NF 


D .TL(L(LNSM ) , SM(LNSV), DII, NFI - 


i TNSV = 1 - NF 


i TKLS = 1 


4. TNSV - INSV + NF 


“^^<LKLS=]KLj 


i - 1 

> calculate L. *m. 

I 11 

; by back substitution 


D 1T..(L(LNSM), SMQNSV), D, NF [ — 


1 = 1 - 

ND = LNSM - 1 



I _ 1 , 

: calculate L 
I by back substitution: 



















E5JL (L, SM, D, NF) 













D(NRW) = (SM (NRW) - S) 






